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Abstract 

Least squares methods based on first-order systems have been recently 
proposed and analyzed for second-order elliptic equations and systems. They 
produce symmetric and positive definite discrete systems by using standard 
finite element spaces, which are not required to satisfy the inf-sup condition. 
In this paper, several domain decomposition algorithms for these first-order 
least squares methods are studied. Some representative overlapping and 
substructuring algorithms are considered in their additive and multiplicative 
variants. The theoretical and numerical results obtained show that the clas- 
sical convergence bounds (on the iteration operator) for standard Galerkin 
discretizations are also valid for least squares methods. 
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1 Introduction 


Least squares methods have been proposed in recent years for second-order el- 
liptic problems, Stokes and Navier-Stokes equations; see Chang [10], Bochev 
and Gunzburger [2], Pehlivanov, Carey, and Lazarov [15], Cai, Lazarov, Man- 
teuffel, and McCormick [5], Cai, Manteuffel, and McCormick [7], Bramble, 
Lazarov, and Pasciak [3], Bramble and Pasciak [4], Carey, Pehlivanov, and 
Vassilevski [8], Cai, Manteuffel, and McCormick [6], Bochev, Cai, Manteuffel, 
and McCormick [1], and the references therein. 

Among the possible approaches, we follow here the one introduced in the 
very recent works of Pehlivanov, Carey, and Lazarov [15] and Cai, Manteuf- 
fel, and McCormick [7]. The second-order elliptic problem is rewritten as a 
first-order system and a least squares functional is introduced. The result- 
ing discrete minimization problem is associated with a bilinear form which is 
continuous and elliptic on an appropriate space. Therefore, the inf-sup condi- 
tion is avoided and standard finite element spaces can be used. The resulting 
linear system is symmetric, positive definite and has condition number of the 
same order as standard Galerkin finite element stiffness matrices, 0(l/h 2 ). 
An interesting alternative approach by Bramble, Lazarov, and Pasciak [3] is 
based on replacing one of the Z 2 -terms in the least squares functional by a 
discrete H~ l - norm. We will not consider here such an alternative. 

The goal of this paper is to extend to these least squares methods some of 
the classical domain decomposition algorithms which have been successfully 
employed for standard Galerkin finite elements. We show that optimal and 
quasi-optimal convergence bounds follow easily from the standard Galerkin 
case. Therefore, domain decomposition provides highly parallel and scalable 
solvers also for first-order system least square discretizations. An overview of 
domain decomposition methods can be found in the review papers by Chan 
and Mathew [9], Dryja, Smith, and Widlund [11], Dryja and Widlund [13], 
and Le Tallec [14]. 

This paper is organized as follows. In the next section, we briefly review 
the first-order system least squares methodology and the main results from 
[7]. In Section 3, we introduce and analyze our domain decomposition algo- 
rithms: overlapping additive Schwarz methods (with coupled and uncoupled 
subspaces; see 3.1), overlapping multiplicative Schwarz methods (3.2), and 
an iterative substructuring method (3.3). In Section 4, we present numerical 
results in the plane that confirm the theoretical bounds obtained. 
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2 First-Order System Least Squares 


We consider the following second-order elliptic problem on a bounded domain 
ft C R 2 or R 3 

' — V • (AVp) + Xp = f in ft, 

< p = 0 on IV, (1) 

n -A Vp = 0 on IV- 

Here A is a symmetric and uniformly positive definite matrix with entries in 
L°°(ft), X is a first-order linear differential operator, IV U Tw = 9ft, and n 
is the outward normal unit vector to IV- 

Defining the new flux variable u = —AVp, the system (1) can be rewritten 
as a first-order system: 


u + AVp = 0 

in ft, 

V • u + Xp = f 

in ft, 

p = 0 

on Tp, 

n u = 0 

on IV, 


This system can be extended to the equivalent system 

u + AVp =0 in ft, 

V • u + Xp = f in ft, 

V x j4 -1 u = 0 in ft, 

p = 0 on Td, 

n • u = 0 on IV, 

k 7 t(A _ 1 u) = 0 on T d , 


( 2 ) 


(3) 


where Vx is the curl operator (in two dimensions Vxu = 0 means ^ — ^ = 
0) and 7 t u = u x n (in two dimensions 7 T u = u • r). 

The following least squares functionals, Go for the system (2) and G for 
the augmented system (3), were studied in [5] ([15] for the case X = 0) and 
[7], respectively: 


Go (v,<7; /) - || v + AVq\\ 2 L 2(Q) + ||V • v + X q — f\\ 2 L 2^ ( 4 ) 

V(v,<?) € Wo (div;Sl) x V, and 

G{v,q‘-f) = H v + + 1| V -v+Xq— /||ii(n) + ||V x (A 1 v)||£ 2 ( fi ) (5) 
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V(v, q) € W x V. 

Here the functional spaces are defined as 


Wo (div, fi) = {v G H(ch'u; 0) : n • v = 0 on T/v}, 

Wo(curlA] f2) = {v G H(cur//1; 0) : 7 t (j4 -1 v) = 0 on To}, 

W = W 0 (div] 0) n W 0 {curlA- : fl), 

V = {q e H\a) : q = 0 on T/j}. 

The least squares minimization problems for (2) and (3) are, respectively: 
Find (u,p) € Wo (div, ft) x V such that 


G' 0 (u,p; /) 


inf 

(v,g)6^Vo(^ l v:^) x V 


G 0 (v,qr;/); 


( 6 ) 


Find (u,p) € W x V such that 


G(u,p;/)= inf G(v j9 ;/). (7) 

(v,g)eWxV 

Simple calculations show that the associated variational problems are, 
respectively: 

Find (u,p) G Wo (div;£l) x V such that 

ao(u, p; v, q) — F(v,q) V(v,g) € W 0 (dit;;ft) x V; ( 8 ) 


Find (u, p) G W x F such that 

a(u,p;v,g) = F(v,g) V(v,g)€ W x V. (9) 

Here the bilinear forms are 


a 0 (u,p; v,g) = (u + 4Vp, v + .AVg)# + (V • u + Xp, V • v + Xq) L 2 , 

<z(u,p;v,g) = a 0 (u,p;v,g) + (Vx (A _ 1 u), V x (A - 1 v)) L a 
and the right-hand side is 

F{v,q) = (/, V- v + Xq) L 2 . 

In [5], it was proved that ao(v, < 7 ; v, < 7 ) is equivalent to (continuous and 
elliptic with respect to) the H(chV, H) x/f^H) norm on Wo(div; Q)xV, under 
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the assumption that a Poincare-Friedrichs inequality holds for p (denoted by 
assumption AO). For the case X = 0, this was already proved in [15]. 

In [7], it was proved that a(v, <7; v, <7) is equivalent to the [H(ch'i;;fl) fl 
H(cur/A;fl)] x H 1 ( fl) norm on W x V, under the same assumption AO. 
Moreover, under three additional technical assumptions denoted by Al, A2, 
A3, it is proven in [7] that a(v,q-v,q ) is equivalent to the H 1 (Q) d+1 norm 
on W x V (d = 2 or 3): 

Theorem 1 Let 6(u,p;v, q) = (u,v)#i( fi )d + (p, q)H l (U) be the bilinear form 
associated with the H 1 ( fl) d+1 norm. 

If the assumptions AO- A3 of [7] are verified, then there exist positive 
constants a and ft such that 

otb(v, q; v, q) < a(v, q-, v, q) V(v, q) € W x V, 

and 

a{u,p:v,q) < 0b{u,p-,u,p) 1/2 b{v,q;v,q) 1/2 
V(u ,p),(v,q) € W x V. 

Because of the equivalence of a(>, •) and b ( •, •), from now on we will concen- 
trate on the variational problem (9) associated with the augmented system 

(3) ' . . ^ 

We introduce a triangulation T/j of fl and associated finite element sub- 
spaces x V h C W x V. We then obtain a finite element discretization 
of (9): 

Find (u^ph) € W^ 1 x V h such that 

a(\ih,ph', Vh,qh) = F(\h,qh) V(vfc,^) € W A x V h . (10) 

For simplicity, we consider continuous piecewise linear finite elements: 

W h = {v <E C°(fl) d : v k \ T € P\(T), VF 6 r^,veW} = Wj 1 x W* x 

v h = {qe C°(Q) : q\ K € P^K), VA' € r h ,q € V). 

and the subscript h for discrete functions will be dropped in the rest of the 
paper. 

Error estimates and results on the conditioning of the resulting stiffness 
matrix can be found in [5] (in [15] for the case X = 0). 

Upon choosing a basis in W k and 14, the discrete problem (10) is turned 
into a linear system of equations Ax = b. We are going to solve such system 
iteratively using domain decomposition techniques. 
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3 Domain Decomposition Algorithms 

We will introduce and analyze our domain decomposition algorithms in the 
Schwarz framework, which has been very successful for standard Galerkin 
finite elements, see [9], [11], [12], [13]. We illustrate the main ideas on algo- 
rithms that are representative of the main classes of domain decomposition 
(additive, multiplicative, overlapping, iterative substructuring). The same 
analysis can be applied to the many other algorithms which have been pro- 
posed and analyzed for the standard scalar case. 

We suppose that the domain is first triangulated by a coarse finite 
element triangulation 77 / consisting of N subdomains of diameter H . The 
fine triangulation 77 * is a refinement of 77 /. For simplicity, we suppose that 
each subdomain is the image under a smooth map of a reference cube. 


3.1 Overlapping Additive Schwarz Methods 

Each subdomain J7, is extended to a larger subdomain $1', consisting of all 
elements of within a distance 6 from $7, (0 < 6 < H). 

Each scalar component of our finite element space W h x V h is decomposed 
as in the standard scalar case: 


N 


N 


N 


N 


w? = £ K 


w 


= £ w h 


2 ,t 


= £ K, V'* = £ V 


rh 


i=i 


i=l 


i=i 


t=i 


where 

= {u € Wj? : support(u) C £)■}, k = 1,2,3, 
V/ 1 = {m£ V h : support (u) C ft-}. 


For each scalar component, a global coarse finite element space is associated 
with the coarse tri angulation 77 / : 


W£ 0 = Wfc — {u € Wf? : u is trilinear on each subdomain f 2*}, k = 1,2,3, 


Vq = V H = {p G V h : p is trilinear on each subdomain 

A first additive method is defined by the following decomposition of the 
discrete space, which maintains the local and coarse coupling between the 
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different scalar components: 


N 

w h X V h = x v ?- 

1=0 


The local spaces are 


Wf x V t h = W* x W*. x Wj,- x Vf i = 1, 2, • • • , N 


and the coarse space is 

w£ x V 0 h = W H xV H = Wf x Wf X wf x V H . 

We define the local projection operators Pi : 'W h x V h —* Wf x Vf by 
a(Pi(u, p); v, q) = a(u, p; v, q) V(v, 9 ) € Wf x Vf, 

and the coarse projection operator To : W ^ 1 x V h — * x Vf by 

a(Po{u,p))\,q) = a(u,p;v,q) V(v,<?) € W£ x V 0 h . 


It is easy to see that the matrix form of the local projections is Pi = 

1 if e* e ft' 

0 otherwise 

trices selecting only the unknowns in ft' for each component and the A{ = 
RiARj are the stiffness matrices of local Dirichlet problems. Analogously, 
To - RhA , where R Jj is the standard piecewise linear interpolation 

matrix from the coarse grid r# to the fine grid for each component, and 
Ah — RhARJi is the coarse grid discretization of our problem (9). Let 


Rj A i 1 jRjA, where the Ri(ek) = 


1 . are the restriction ma- 


tt 

Paddl = Y: Pi- 

i = 0 

The original discrete problem is then equivalent to the preconditioned prob- 
lem 

Paddl (u, p) = Qaddl 

where g = £iIo T(u,p) ; see Chan and Mathew [9]. In matrix form, this 
problem can be written as M~ l Ax = M~ l b , where the preconditioner is 
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M- 1 =Ef =1 Rj A~ l Ri + RJiAJj 1 Rtf. An optimal convergence bound for this 
algorithm is given in Theorem 2. 

A second additive method is obtained by dropping the coupling between 
the different scalar components of u and p. Uncoupled local spaces are now 
defined bv 

W‘, = n;‘, x {0} x {0} x {0}, 
w‘ f = {0} x W* x {0} x {0}, 
w$, = {0} x {0} x Wl x {0}, 

v* = {o) x {o} x {o} x v;\ 

and the coarse spaces by 

w" = W‘ 0 = VK» 0 x {0} x {0} x {0}, 

W" = W*„ = {0} x x {0} x {0}, 

W? = W* 0 = {0} x {0} x W£„ x {0}, 

\ H = v£ = {0} x {0} x {0} x V 0 h . 

We then have the following decomposition 


N 


N 


N 


N 


w» x = X Wf,, + E w t, + E w ‘, + E v .‘ + w ? + W 2 H + W 3 H + v 


H 


2 = 1 


1 = 1 


1=1 

3 N 


i=l 


EE w t + E v t 


k= 1 i=0 


N 

C 

1=0 


As before, we define projections Pk,i : W* x V h — *• W Ar = 1,2,3, i = 
0, !,•••, N and x V h — > Vf, i = 0, N. and the additive 


operator 


3 N N 

Padd2 = E E a,, + E p*.<- 


fc = 1 1=0 1=0 

We note that this algorithm can equivalently be defined by the same choice 
of subspaces as for P a dd\ but using the bilinear form &(•,•) (introduced in 
Theorem 1) instead of g(-,-) in the definition of the projections. In fact this 
uncoupled preconditioner corresponds to applying four identical copies of a 
scalar preconditioner to each scalar component. An optimal bound also holds 
for this algorithm. 



Theorem 2 There exists a positive constant C independent of h,,H and 6 
such that 

cond(P) < C{ 1 + ^), 

0 

where P = P a ddi or P — P add2 . 

Proof. An upper bound on the spectrum of P is standard, since each point of 
fl belongs to a fixed number of extended subdomains independent of N (for 
example, for 6 < H/2 each point belongs to at most four (in 2D) or eight (in 
3D) extended subdomains). A lower bound is obtained by classical Schwarz 
analysis. 

For P = Paddi , since we use exact projections, the theorem is equivalent 
to the following partition property (see Dryja and Widlund [13] or Chan and 
Mathew [9]): 

There exists a constant Co such that V(u,p) € x V h , there exists a 
decomposition (u,p) = X^ 0 (u t ,p;), w ith (u ,,p,) € W f x Vf such that 

N 

5^a(u„p t ;u„pi) < C 0 2 a(u,p;u,p). 

1=0 

By the equivalence of Theorem 1, this inequality is equivalent to 
N 

X! l( u iiP«)lrn»w-H ^ C’ 0 2 |( u ,p)|J H i ) d + i, 

i=0 

which is a direct consequence of the scalar result proven by Dryja and Wid- 
lund [13]: 

XI — Cl\u k \ 2 H i, XI IPilw 1 ^ CoIpIh 1 ? 

t=0 t=0 

with Cl = C{ 1 + f ). 

For P = P a dd 2 1 since the subspaces are the same but we use inexact 
projections defined by b ( •,•) instead of a(-,-), we need only to show that 
there exists a constant u such that a(u,p;u,p) < a;6(u,p;u,p) V(u,p) € 
Wf x Vf, i — 0,1 (see Dryja and Widlund [12]). This follows 

immediately from the equivalence of a and b. 
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3.2 Overlapping Multiplicative Schwarz Methods 

By using the same coupled local and coarse spaces as in the additive algorithm 
Padd\i we can define a multiplicative operator: 

Pmult =/-(/- P N ) •••(/- />,)(/ - P 0 ). 

The multiplicative algorithm consists in solving the nonsymmetric system 

Pmultiy^iV) — S mult 

by an iterative method such as GMRES. 

We can also define a symmetrized multiplicative operator 


Pmults = /-(/- Po) •••(/ - Pn-i)(I ~ Pn)(I ~ P/V-l) •••(/ — P 0 ) 


and a symmetrized algorithm, consisting in solving the symmetric system 

P multsi^-i P/ — pmults 

by an iterative method like CG. 

Theorem 3 There exists a positive constant C independent of h, H and 6 
such that 

cond(P mu i ts ) < C{\ + —). 

The proof is again based on the extension of the scalar result (see, for exam- 
ple, Chan and Mathew [9]) by using the equivalence of Theorem 1. Analo- 
gously, multiplicative versions of P a dd2 could be built using uncoupled local 
and coarse spaces. 

3.3 An Iterative Substructuring Method 

For a complete and detailed analysis of this class of methods, we refer to 
Dryja, Smith and Widlund [11]. Here we consider only a simple represen- 
tative of this class, namely the analog of Algorithm 6.2 in [11], which is 
vertex-based and has a standard coarse space. For simplicity, we only con- 
sider the uncoupled additive version. 
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The standard first step of nonoverlapping methods is the elimination of 
the variables interior to each subdomain (at least implicitly). We then work 
with the Schur complement S = Kbb - I<J b K ii K ib of the stiffness matrix 


I\ = 


Kii Kib \ 
KJb Kbb ) ' 


The reduced linear system with S involves only variables on the interface 
T = Udfl,- \ Tp. When solving with a preconditioned iterative method, we 
need only the action of S' on a given vector and there is no need to assemble 
S explicitly. 

In the Schwarz framework, working with S corresponds to working with 
the discrete harmonic subspace W h x V h of the original space W /l x V h . 
Local spaces are associated with the geometric objects (faces F t . edges E t 
and vertices t?,) forming the interface T. Each scalar space is decomposed as 

K = E Kf, + E Kb, + E K„, k = 1 . 2 , 3 , 

F x E x v i 

and 

v"‘ = E ' , 4 + E^+EKJ- 

F t Ei v \ 

Here, for example, W£ F . = {u € W h : u = 0 on T* — F,^}, where T/, and F Xi h 
are the set of nodes on T and F, respectively. The other spaces are defined 
analogously. As for the overlapping case, we then embed these scalar spaces 
in our product space ~W h x V h : for example, Wf x {0} x {0} x {0}- 

As a coarse space, we consider the discrete harmonic subspace of the same 
coarse space used for P a dd 2 , i.e., + + Y H . We obtain the 

following decomposition 


W l x v k = E(E wjU + E K,e, + EwE,„ + wf > 


*=1 Fi 


E t 


+ E Vf. + E v|, + E vi + v". 

F, E, «. 

By defining as before projection operators into the subspaces, we form the 
additive operator 

Pi, = E(E a.f. + E ru + E p*.«, + Pk o). 

k=\ F, E, v t 


10 



where again for k — 4 the projections are into the V* spaces. 

Theorem 4 There exists a positive constant C independent of h and H such 
that 

cond(P ls )<C( l+\og(H/h)) 2 . 

As before, the proof is based on the extension of the scalar result (see Dryja, 
Smith and Widlund [11], Theorem 6.2) by using the equivalence of Theorem 
1. 

4 Numerical Results 

In this section, we report the results of numerical experiments which confirm 
the optimal convergence bounds obtained in the previous sections. All the 
results have been obtained with Matlab 4.2 running on Sun Sparcstations. 
The model problem considered is the standard Poisson equation ( A — I, X = 
0) on the unit square, with p = 0 on To = dTl and 7 r u = 0 on dQ (i.e. U\ = 0 
on {y = 0} and {y = 1}; u 2 = 0 on {x = 0} and {x = 1}). The right-hand 
side / is chosen such that we have p(x, y) = sin(7rx)sin(iry) as exact solution, 
fi is decomposed into a regular grid of N square subdomains, with N varying 
from 2 x 2 to 8 x 8. The fine grid mesh size h varies from 1/32 to 1/128. 

The Krylov method used for all the symmetric problems is PCG, while 
we use GMRES for the nonsymmetric problem with P mu it ■ The initial guess 
is always zero and the stopping criterion is | || 2 / j]x-o H2 < 10 -6 , where 77. is 
the residual at step k. 

The local and coarse problems involved in the application of the precon- 
ditioners are always solved directly. For each method, we report the number 
of iterations and Lanczos-based estimates of the condition number and the 
extreme eigenvalues (except for the multiplicative algorithm, where we report 
the average convergence factor instead). 

Overlapping additive methods. We have first studied the coupled method 
Paddi with fixed minimal overlap size 6 = h. The mesh size h is decreased 
while the number of subdomains N is increased proportionally, so that the 
subdomain size H/h = 16 is kept constant ( H — l/y/N). The results are re- 
ported in Table 1 and clearly show a constant condition number cond(P a< idi ) = 
A ma x/A m in, for problem sizes from 3007(A r = 4) to 48895(A r = 64). 
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Table 1: P a ddi ■ Overlapping Additive Schwarz with fixed overlap size 6 = h. 


N 

ESI 

iter. 

cond( Paddi) 

^ max 

^ min 

4 

32 

16 

11.2172 

4.0048 

0.3570 

9 

48 

19 

12.1787 

4.0068 

0.3290 

16 

64 

20 

11.9775 

4.0050 

0.3343 

25 

80 

20 

11.1689 

4.0052 

0.3586 

36 

96 

21 

12.5450 

4.0044 

0.3192 . 

49 

112 

20 

11.9944 

4.0050 

0.3339 

64 

128 

21 

12.5500 

4.0047 

0.3191 


Table 2: P a ddi • Overlapping Additive Schwarz with fixed number of subdo- 
mains N — 64. 


8 

h 1 

iter . 

cond(P ad di) 

^max 

^ min 

h 

128 

21 

12.5500 

4.0047 

0.3191 

2h 

128 

17 

7.1316 

4.0307 

0.5651 

3h 

128 

16 

5.5769 

4.0765 

0.7309 

4h 

128 

15 

4.9540 

4.1396 

0.8356 

5h 

128 

15 

4.6460 

4.2170 

0.9076 

6h 

128 

15 

4.5125 

4.3054 

0.9541 

7h 

128 

16 

4.5859 

4.4018 

0.9598 
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Table 3: P a dd2 '■ Overlapping Additive Schwarz with fixed overlap size 8 = h. 


N 

h' 1 

iter . 

cond(P add2 ) 

^max 

^ min 

4 

32 

17 

10.3521 

4.0050 

0.3868 

9 

48 

20 

12.6290 

4.0051 

0.3171 

16 

64 

20 

11.9811 

4.0051 

0.3342 

25 

80 

21 

11.3821 

4.0052 

0.3518 

36 

96 

21 

12.5458 

4.0043 

0.3191 

49 

112 

20 

11.9997 

4.0052 

0.3337 

64 

128 

21 

12.5261 

4.0047 

0.3197 


Table 4: P a dd2 ■ Overlapping Additive Schwarz with fixed number of subdo- 
mains N = 64. 


r u 

hr 1 

iter . 

cond(P add2 ) 

max 

min 

"h 

128 

21 

12.5261 

4.0047 

0.3197 

2h 

128 

17 

7.1206 

4.0315 

0.5661 

3h 

128 

16 

5.5513 

4.0777 

0.7345 

4h 

128 

16 

5.3850 

4.1442 

0.7695 

5h 

128 

16 

5.4545 

4.2233 

0.7742 

6h 

128 

16 

5.5306 

4.3158 

0.7803 

7h 

128 

16 

5.6176 

4.4297 

0.7885 
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Table 5: P m ult and P mu its : Overlapping Multiplicative Schwarz with fixed 
overlap size S = h. 


N 

i k- 1 

multiplicative (GMRES) 
iter. p = (n/ro) 1 ^ 

symmetrized multiplicative (CG) 
iter. cond(P mu i ts ) A max A min 

4 

32 

8 

0.1847 

7 

1.8576 

0.9994 

0.5379 

9 

48 

7 

0.1433 

6 

1.7398 

0.9999 

0.5749 

16 

64 

6 

0.1233 

6 

1.7600 

0.9999 

0.5681 

25 

80 

6 

0.1102 

6 

1.6810 

0.9999 

0.5948 

36 

96 

6 

0.1021 

6 

1.6940 

0.9999 

0.5902 

49 

112 

6 

0.0952 

6 

1.6661 

0.9999 

0.6001 

64 

128 

5 

0.0849 

6 

1.7308 

0.9999 

0.6079 


In Table 2, we fix the mesh size ( h = 1/128) and the decomposition 
( N — 64) and we vary the overlap size 6 from h to Ih. As in the scalar case, 
the condition number cond{P a dd\ ) improves as 6 increases, because of A m ; n 
being closer to unity. For large overlap, the improvement becomes negligible 
or negative, because of the growth of A max . 

The same sets of results for the uncoupled method P a d<n are reported 
in Table 3 and Table 4, respectively. For this simple model problem, the 
uncoupled method is only slightly worse than the coupled one, in terms of 
iteration count (some condition numbers are almost the same or even better 
for P a dd 2 ) • We point out that although A = /, eliminating diffusive coupling 
between the flux components, there is still coupling between the flux variables 
and p, so the strong performance of P a dd 2 is encouraging. 

Overlapping multiplicative methods. In Table 5, we compare the multi- 
plicative method P mu it accelerated with GMRES and the symmetrized mul- 
tiplicative method P mu it s accelerated with CG. We consider the two methods 
with minimal overlap and constant subdomain size. Since P mu it is nonsym- 
metric, we report the average convergence factor p = {r t jr G ) 1 ^ instead of the 
condition number. Even if the symmetrized version is approximately twice 
as expensive as the standard one, the number of iterations is almost the same 
for the two methods. Therefore, the symmetrized version is less efficient on 
this simple problem. 
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Table 6: Pj£ Iterative Substructuring . 


N 

h - 1 

iter. 

cond(P ls ) 

^max 

min 

4 

32 

9 

3.4035 

1.5691 

0.4610 

9 

48 

17 

7.8812 

1.8497 

0.2347 

16 

64 

18 

7.8543 

1.7962 

0.2287 

25 

80 

18 

8.5822 

1.8864 

0.2198 

36 

96 

19 

9.4115 

1.8511 

0.1966 

49 

112 

18 

8.6646 

1.8939 

0.2185 

64 

128 

19 

9.6532 

1.8617 

0.1928 


Iterative sub structuring. Table 6 shows the results for the iterative sub- 
structuring methods P zs with fixed subdomain size. They clearly show a 
constant bound for the condition number and the number of iterations. 


5 Conclusions 

In this paper, some domain decomposition algorithms have been introduced 
for the discrete systems arising from first-order system least squares methods 
applied to second-order elliptic problems. These recently proposed methods 
allow the use of standard finite element spaces, which are not required to 
satisfy the inf-sup condition. 

The analysis of the domain decomposition algorithms follows from analo- 
gous results for the standard Galerkin case and the equivalence between the 
bilinear form associated with the least squares functional and the H 1 ( fi) d+1 
norm. 

Optimal convergence bounds have been proven for overlapping algorithms 
(additive, multiplicative, coupled, uncoupled versions), while quasi-optimal 
bounds have been proven for iterative substructuring algorithms. Numerical 
experiments on a simple model problem confirm these bounds. 

Future work will investigate the performance of these algorithms for prob- 
lems with convection and for elliptic systems. 
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